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Abstract 

The two dimensional Hubbard model with a single spin-up electron interact- 
ing with a finite density of spin-down electrons is studied using the quantum 
Monte Carlo technique, a new conjugate gradient method for the evaluation 
of the Edwards wavefunction ansatz, and the standard second order pertur- 
bation theory. We performed simulations up to 242 sites at U/t = 4 reaching 
the zero temperature properties with no "fermion sign problem" and found a 
surprisingly good accuracy of the Edwards wavefunction ansatz at low density 
or low doping. The conjugate gradient method was then applied to system 
up to 1922 sites and infinite U for the Edwards state. Fermi liquid theory 
seems to remain stable in 2D for all cases studied with the exception of the 
half filling case where a "Luttinger like behavior" survives in the Hubbard 
model , yielding a vanishing quasiparticle weight in the thermodynamic limit. 
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The anomalous properties found in several High T-c oxides have renewed a considerable 
attention and an increasing interest for the simplest models which may explain non con- 
ventional behavior beyond the Fermi-liquid theory of normal metals and the BCS theory 
of usual superconductors.. The Hubbard model is now the most popular hamiltonian in 
condensed matter physics and the search for a satisfactory solution in two spatial dimension 
represents one of the most important challenge in theoretical physics. In one dimension we 
now have a complete solution, not only for the energy spectrum, but also for asymptotic 
properties of correlation functions. The physics of the ID Hubbard model is well described 
by the Luttinger liquid theory. The Fermi liquid theory is unstable in this case due to di- 
vergences in perturbation theory (PT). Such divergences are usually not present in higher 
dimensionality and a possible anomalous phase in 2D - as suggested by P. W. Anderson |J 
- can be explained only within a non perturbative approach. 

Historically before the Lieb-Wu exact solution || a much simpler case was solved ; it is 
the case when only one particle with spin up interacts with a finite density of spin down 
electrons. Although this problem is very much simplified and probably far from reality, 
it already contains the basic features of ID conductors : the quasiparticle weight vanishes 
in the thermodynamic limit, and the spectrum consists of holon and spinon elementary 
excitations ||£§ 

In this paper we attempt to search for a non Fermi liquid phase in the 2D Hubbard model 
in this simplified sector. In fact although we do not provide an exact analytical solution, 
the numerical advantages will be clearly evident in this single spin- flip Hilbert-space of the 
model. 

We thus consider the Hubbard model on a bipartite lattice , where the spin-up and a 
finite density p = j- of iV spin down electrons hop trough the L lattice sites with different 
hopping amplitudes, t c and t d respectively. The hamiltonian therefore reads: 

H = X) *c 4+r M CR + U 4+r M d R + U X C R°R d R d R ( 1 ) 

R,t„ R 

where d (c) are fermion operators for spin down (spin up) electrons, r M are the 2d nearest 
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neighbor vectors and d (= 2 ) is the spatial dimension. 
We are interested in the quasiparticle weight : 

Z P = \ < i> G \dl\il>F > I 2 - (2) 

where > is the ground state of H in the single spin flip subspace and with finite momen- 
tum p, and \ipF > = II c l|0>is the non interacting electron sea of spin down electrons, 
i.e. the ground state of H without d— electrons. Such a state is assumed in the following 
to be a non degenerate , translation invariant closed shell state where = t c J2 e ,and 
ep is the Fermi energy of the down-spin electrons, [ipp > satisfies these properties only for 
particular number of electrons, e.g. N = 1, 5, 13 . . . ., covering all possible densities in the 
thermodynamic limit. 

As well known, the quasiparticle weight Z p measures the strength of the 5— function in 
the spectral function and is finite in the thermodynamic limit if Landau- Fermi liquid theory 
is valid. For example in one dimension, where Fermi liquid theory breaks down, Z vanishes 
as a power law in the thermodynamic limit ,as predicted by the 1 Luttinger liquid theory 
Z oc N~ e . The Hubbard model in the single spin-flip sector is consistent with this general 
solution |3j and thus represents one of the simplest toy models where Fermi-liquid theory 
can be tested in higher dimensions. 

The hamiltonian H can be further simplified by tracing out exactly the d— electron com- 
ponent from the ground state wavefunction. In fact, using translation invariance we can 
write the ground state \ipQ > of total momentum p in the following form: 

\^g>= 4fE e " ip * T ^ohA> (3) 

Viv R 

where \ip > is a wave function depending only on the spin- down electrons, Tr is the 
translation operator by a vector R (T R Cr> T_ r = c R+R ') and O is the origin of the lattice. 
Then the following effective hamiltonian H' for \ip > is obtained: 

H' = t c J2 Cr+^cr + T Tm + U cW (4) 
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Note that H' is not translation invariant and that the on site Coulomb repulsion U becomes 
now a simple one-body contribution, i .e. quadratic in the fermion fields. Using Eq.([3|) Z p 
in Eq.(Q) is replaced by the overlap between the U = {\ipF > )and the finite U ground 
state, at fixed number of electrons : Z p = | < ipF^ > \ 2 ■ Thus the question of Fermi liquid 
or non Fermi liquid theory in this model is simply related to the stability of the ground state 
of H' under the local perturbation U CqCo- 

Another useful quantity which we will consider in the following is the momentum distri- 
bution of the d— electron: rik = < ipG\d\dk\ipc > ■ By means of Eq. fl3|) n& is related to the 
expectation value of the p — fc-momentum projector on the state \ip >: rik = < ip\P p _k\4' >> 
where the projector on the subspace of momentum Q is Pn = jJ2 e~ l< ^ R T^. By inserting 

R 

in < i[)\Pp^k\i/; >, a complete set > of translation invariant states containing >, it 
easily follows that: 

n p > Z p (5) 

The previous relation can be viewed as a particular case of the Migdal theorem relating 
the jump of the momentum distribution at the Fermi surface to the amplitude of the spec- 
tral weight. In fact we expect that the inequality (|5]) turns in an exact equality in the 
thermodynamic limit. 

As it is easy to verify the ground state of the hamiltonian (|4]) is a free electron Slater 
determinant in several limiting cases. For N = 1 -corresponding to the two electron problem 
for the hamiltonian H- there is of course no correlation in H'. For t^ = the hamiltonian 
becomes the well known Falikov-Kimball model and the effective hamiltonian H' is free and 
exactly solvable. The ground state in presence of the local perturbation Uc^cq is orthogonal 
to the non interacting state, yielding Z p cx N~ e 0. Finally for U = as well as for tj — > oo 
the free Fermi gas tpp is the ground state and Z p — 1. 

The above limiting cases are not surprising since only the term proportional to t^-the 
spin-up kinetic term-is a true many body term in the effective hamiltonian H', all remaining 
ones being one-body contributions. Moreover this spin up kinetic term is obviously not 
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extensive in the size of the system and can weakly affect the correlation in the ground state. 
Based on the previous considerations it is likely that the ground state of H' is always very 
close to a simple Slater determinant and thus an Hartree-Fock (HF) wavefunction (i.e. the 
Slater determinant which minimize the expectation value for the energy) may have a very 
good overlap with the exact ground state of H' . 

The HF wavefunction of the hamiltonian (1) \ipEWA > is nothing but the Edwards 
wavefunction ansatz (EWA), which is exact in ID for td = t and very accurate in energy 
in the 2D case ||. Nevertheless such an Hartree-Fock wavefunction corresponds ,by Eq.(f|) 
, to a non trivial correlated state of the Hubbard hamiltonian H. It is actually a linear 
combination of L- free electron states, yielding for instance the Bethe -ansatz wavefunction 
in ID 

We have used a well known Quantum Monte Carlo (QMC) technique [[7] to evaluate the 
quasiparticle weight Z p and the momentum distribution n p in this simple model for p = 
and t c = td- ■ The ground state \ipc > is filtered out by imaginary time propagation of 
a given trial wavefunction \%l) T >, \ip G >cx lim e~ tH \ipT >, after the usual Trotter-Suzuki 
decomposition of the imaginary time in P slices of length At = 4-|'0T >= l^F > i n a U 
the present simulations. The unrestricted Hartree-Fock calculations for the EWA and the 
straightforward second order perturbation theory for Z p and n p -coinciding at this order- 
, are compared with the QMC simulations on finite lattices , ly/2 x ly/2, with periodic 
boundary conditions tilted by 45 degrees and odd I. The convergence in imaginary time is 
systematically reached within statistical errors in all cases studied. In fact for the Hubbard 
model in the single spin-flip sector the QMC does not suffer the so called "fermion sign 
problem", since in the worst case (I = 11) the average sign is approximately 0.6. 

In all the QMC simulations I have always found n p=0 and Z p=0 equals within statistical 
errors, with an error bar for n p three times smaller than the one for Z p . n p is an upper 
bound for the quasiparticle weight (Eq.|^) but it is always very close to Z p . For instance 
when td = or within the EWA approximation n — Z < 10~ 4 for all sizes studied. In the 
following we thus identify the two quantities for the sake of simplicity. 
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We got very accurate results for the evaluation of the EWA wavefunction using a new 
conjugate gradients (CG) technique for electronic structure calculations. In order to apply 
the CG algorithm H to the minimization of an energy functional E({ip n }) = ^^~!^r^j~^~ > 
the orbitals ip n {k) of the Slater determinant \ip > have to be orthogonalized from time 
to time for numerical stability ||. In doing so one spoils the efficiency of the conjugate 
directions, and the minimization of the energy becomes slower. This is why the standard 
CG algorithm with orthogonalization does not improve much the steepest descent scheme, 
as also discussed in [|J. 

In order to solve the previous difficulty we have used the following simple scheme. Since 
the functional E({ip n }) is invariant for any transformation of the orbitals ip n — > Y^m A. n , m i\) m 
it is possible to choose the NxN matrix A such that the orbitals read: ip n (k) = 6k,k n +u n (k), 
where k n are wavevectors inside the Fermi surface < ep and the functions u n (k) all 
vanish for tk < e^. The number of independent degrees of freedom is thus N(L — N), as it 
should be from general ground. In this way we fix the gauge of the transformations, that 
leave E unchanged, and we apply the CG strategy, without orthogonalization requirements. 
The method is of course not restricted to the plane wave basis for the orbitals |J. In the 
largest sizes simulations (~ 2000 sites) this novel CG minimization is approximately an order 
of magnitude more efficient compared to the usual steepest descend method. This factor 
increases with the size of the system, opening new possibilities for large scale simulations. 
Further details of this new method will be published elsewhere ||. 

At low density we have studied closed shell systems with the closest density just below 
and above the value p = \- This sequence should converge to p = j in the thermodynamic 
limit and should minimize size effects. In Fig. p] we see that size effects are very important 
and it is not possible to obtain reasonable conclusions with a small size calculation. The 
EWA is exact for N — 1 (smallest size in Fig. [I] ) and is a test of our QMC scheme. We note 
that the results obtained with the EWA are practically indistinguishable from the QMC 
data, yielding a strong support for the accuracy of the EWA even for relatively large size 
(up to 242 sites), not accessible by exact diagonalization. 



Further evidence for the accuracy of the EWA is given by the explicit calculation of the 
energy fluctuations AE 2 =< ijj\H 2 \ijj > - < ip\H\vp > 2 on the EWA \tp >= \ip E wA > 
AE 2 vanishes for an exact eigenstate and is given by U 2 p(l — p) for the Fermi gas wave- 
function \ifjp >• The Edwards wavefunction typically improves this variance by three order 
of magnitude compared to the Fermi gas wavefunction which still is very good for this small 
positive U value. As shown in Fig. [I] the behavior of the QMC data up to I — 11 and the 
EWA data for larger systems up to I = 31 are very similar to the PT, which is finite for 
d = 2 and p ^ |. . Furthermore EWA and QMC data are always well above the PT results 
, strongly suggesting that Z p should be finite in the exact calculation as well as for larger U 
{U — oo is shown in Fig. [I] ) . Similar scenario appears also evident for negative U (Fig. 0) 
, corresponding by the particle -hole transformation to positive U at density 1 — p = | , 
i.e. the low doping region. However in this case the PT results for Z p are larger than the 
QMC and EWA data and the conclusion of a finite Z p , although quite likely, is less clear. 
Note also that for infinite U, Z p seems to drop at the largest sizes. 

At half filling L/2 = N = I 2 the U = ground state is a non-degenerate closed shell. The 
Fermi surface is a square and is commensurate with the finite mesh in the Brillouin zone, 
leading to size effects very smooth and well controlled. In Fig.[3] we show Quantum Monte 
Carlo results as a function of At 2 . The At 2 — > limit should be considered a formally exact 
and unbiased property of the ground state for long enough t. The results obtained with the 
EWA are systematically larger than the QMC data, except for small sizes. 

In this case the second order PT is logarithmically divergent due to the nesting of the 
Fermi surface: \nZ p=0 = -f/ 2 0.00335002 ln(L) + 0(U 4 ). As in ID this kind of divergence 
suggests a power law decaying for Z p oc N~ e . In fact as shown in Fig.[3](b) the ln-ln plot 
appears linear for ln(L) — ► oo both for the PT data {Opt = 0.0536002) and for the EWA 
ones {9 ewa — 0.028). The QMC data lie in between the two straight lines giving a strong 
evidence that Z p eventually vanishes as a power law wit 6 ~ 0.04. This is maybe one of the 
first example of Luttinger-like behavior in a 2D system and represents the central result of 
this paper. It is also interesting to note that the Luttinger liquid exponent for U = oo is very 
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close to 9 = | for the EWA data, exponent well known in ID where it can be determined 
exactly [p~0| ,p7T|] . 

In conclusion we have presented an accurate and well controlled size scaling of the quasi- 
particle weight in a numerically tractable sector of the 2D Hubbard model. We have used 
QMC and a new conjugate gradient algorithm for the evaluation of the EWA. This ansatz 
turns out to be an almost exact approximation in most cases for U = \4t\. 

The quasiparticle weight Z p looks always finite with the noticeable exception of p = |, 
where the nesting of the Fermi surface leads to a power law decaying Z p , as well as to a 
logarithmically divergent PT. At low doping and large U a non perturbative break-down of 
Fermi liquid theory is not inconsistent with our data , and further study is needed to clarify 
this issue. 

I acknowledge useful discussions with A. Parola, M. Fabrizio and E. Tosatti. Use of the 
Cray-YMP was supported by CNR, project "Sistemi Informatici e Calcolo Parallelo" . 
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FIGURES 

FIG. 1. (a) QMC data (empty dots) vs. EWA ones (full squares) for U = 4i c and t c = td- The 
continuous line (upper p < 1/4) connects EWA data for N = 1,9,21,37, the lower curve (dashed 
p > 1/4 ) for N = 5, 13, 25, 45, 61. The QMC data are converged in imaginary time for t = I + 4 
with / = 3, 5, 7, 9, 11 from left to right. The At correction was estimated for / = 3, 5, 7 with several 
points (see Fig. ^). (b) EWA data for U = 4i c (continuous line) and U = oo ( long dashed line) up 
to Z = 31 compared with second order PT results (dashed line). The arrow indicates the infinite 
size PT result. The lines are guides to the eye. 

FIG. 2. same as Fig.l for negative U. The QMC data refers to I = 3, 5, 7, 9. Similar consider- 
ations apply for the At corrections and the imaginary time error was negligible for t = I + 6. 

FIG. 3. (a) logarithmic plot of the quasiparticle weight. The empty dots are QMC data after 
extrapolation to At — ► 0. The continuous line and the long dashed line connect EWA data for 
U = 4t c ( full dots) and U = oo (full squares) respectively. The number of electrons was fixed to 
N = I 2 for I = 3, 5, . . . 23. The PT data (full triangles) were calculated up to I = 53 (not shown in 
the picture) and the dashed line is the exact slope in the thermodynamic limit, (b) QMC data as 
a function of At 2 for U = 4i c and the imaginary time t = 1 + 6, large enough to have converged 
results. The continuous lines are least square fit of the data and the dashed line (I = 9) has a slope 
estimated from the smaller sizes. The arrows indicate the EWA values for I = 3, 5, 7, 9 from top to 
bottom figure respectively. 
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